
#ifndef _GPUSORT_KERNEL_H_
#define _GPUSORT_KERNEL_H_

#include <stdio.h>
#define TPB 512

__device__ inline void swap(float & a, float & b)
{
	// Alternative swap doesn't use a temporary register:
	// a ^= b;
	// b ^= a;
	// a ^= b;

    float tmp = a;
    a = b;
    b = tmp;
}
__global__ static void mergesort(float * values, int dimx, int dimy, int dimz)
{
	__shared__ float sharedA[TPB];

	const unsigned int tid = threadIdx.x;
	const unsigned int bidx = blockIdx.x;
	const unsigned int bidy = blockIdx.y;

	int blocksPerColumn=(dimy/TPB);

	const unsigned int idk= bidy;
	const unsigned int idi= bidx/blocksPerColumn;
	const unsigned int idj= (bidx - idi*blocksPerColumn) * TPB + tid;


	// Copy input to shared mem.
	sharedA[tid] = values[idk*dimx*dimy + idj*dimx +idi];

	__syncthreads();

	// Parallel bitonic sort.
	for (unsigned int k = 2; k <= TPB; k *= 2)
	{
		// Bitonic merge:
		for (unsigned int j = k / 2; j>0; j /= 2)
		{
			unsigned int ixj = tid ^ j;

			if (ixj > tid)
			{
				if ((tid & k) == 0)
				{
					if (sharedA[tid] > sharedA[ixj])
					{
						swap(sharedA[tid], sharedA[ixj]);
					}
				}
				else
				{
					if (sharedA[tid] < sharedA[ixj])
					{
						swap(sharedA[tid], sharedA[ixj]);
					}
				}
			}

			__syncthreads();
		}
	}

    // Write result.
    values[idk*dimx*dimy + idj*dimx +idi] = sharedA[tid];

}

__global__ static void mergeblocks(float * values, float * tvalues, int k, int dimx, int dimy, int dimz)
{
	const unsigned int tid = threadIdx.x;
	const unsigned int bidx = blockIdx.x;
	const unsigned int bidy = blockIdx.y;

	const unsigned int blocksPerColumn=dimy/TPB;

	const unsigned int idk= bidy;
	const unsigned int idi= bidx/blocksPerColumn;
	const unsigned int idj= (bidx-idi*blocksPerColumn)*TPB + tid;

	int lmid = ((int)((bidx-idi*blocksPerColumn)/k)*2 +1)*TPB *k/2;
	int ltid= idj;

	int gtid= idk*dimx*dimy + ltid*dimx +idi;

	int newValue = values[gtid];
	int newltid;

	__syncthreads();

	if(ltid<lmid)
	{
		int pos=0;
		int rankA= ltid;

		//We need to do binary search in the second half of the block
		int tidStart=lmid;
		int tidEnd= lmid + (k/2)*TPB -1;

		//Binary search from tidStart to tidEnd to find value

		int low= tidStart;
		int high=tidEnd;
		while (low<=high)
		{
			int lmidb= (high+low)/2;
			int midValue=values[idk*dimx*dimy + lmidb*dimx + idi];

			if(newValue<midValue)
			{
				high=lmidb-1;
			}
			else if (newValue>midValue)
			{
				low= lmidb+1;
			}
			else if (newValue==midValue)
			{
				while (values[idk*dimx*dimy + lmidb*dimx + idi]== newValue && lmidb>=tidStart)
				{
					lmidb--;
				}
				lmidb++;
				pos=lmidb-lmid;
				break;
			}

			pos=low-lmid;
		}

		newltid= rankA+pos;
	}
	else
	{
		int pos=0;
		int rankB= ltid- lmid;

		//We need to do binary search in the first half of the block
		int tidStart=lmid- (k/2)*TPB;
		int tidEnd=lmid-1;

		//Binary search from tidStart to tidEnd to find value

		int low= tidStart;
		int high=tidEnd;
		while (low<=high)
		{
			int lmida= (high+low)/2;
			int midValue=values[idk*dimx*dimy + lmida*dimx +idi];

			if(newValue<midValue)
			{
				high=lmida-1;
			}
			else if (newValue>midValue)
			{
				low= lmida+1;
			}
			else if (newValue==midValue)
			{
				while (values[idk*dimx*dimy + lmida*dimx + idi]== newValue && lmida>=tidStart)
				{
					lmida--;
				}
				lmida++;
				pos=lmida;
				break;
			}

			pos=low;
		}

		newltid= rankB+pos;
	}

	tvalues[idk*dimx*dimy + newltid*dimx +idi] = newValue;

}


__global__ static void swapvectors(float* values, float* tvalues, int dimx, int dimy, int dimz)
{
	const unsigned int tid = threadIdx.x;
	const unsigned int bidx = blockIdx.x;
	const unsigned int bidy = blockIdx.y;

	const unsigned int blocksPerColumn=dimy/TPB;

	const unsigned int idk= bidy;
	const unsigned int idi= bidx/blocksPerColumn;
	const unsigned int idj= (bidx-idi*blocksPerColumn)*TPB + tid;
	values[idk*dimx*dimy + idj*dimx +idi]= tvalues[idk*dimx*dimy + idj*dimx +idi];
	__syncthreads();
}

__global__ static void cloneholes(float* values, float* tvalues,int k, int dimx, int dimy, int dimz)
{

	const unsigned int tid = threadIdx.x;
	const unsigned int bidx = blockIdx.x;
	const unsigned int bidy = blockIdx.y;

	const unsigned int blocksPerColumn=dimy/TPB;

	const unsigned int idk= bidy;
	const unsigned int idi= bidx/blocksPerColumn;
	const unsigned int idj= (bidx-idi*blocksPerColumn)*TPB + tid;

	int tempidj=idj;

	while(values[idk*dimx*dimy + tempidj*dimx +idi]==-1)
	{
		tempidj--;
		if(values[idk*dimx*dimy + tempidj*dimx +idi]!=-1)
		{
			values[idk*dimx*dimy + idj*dimx +idi]=values[idk*dimx*dimy + tempidj*dimx +idi];
			break;
		}
	}


	tvalues[idk*dimx*dimy + idj*dimx +idi]= -1;
	__syncthreads();

}
#endif // _GPUSORT_KERNEL_H_
